Projet

##Examen des Données:

Avant de se plonger dans des analyses spécifiques, il est essentiel d’examiner les données pour comprendre leur structure, leurs variables et leurs statistiques de base. Nous chargerons les données, examinerons les premières lignes et passerons en revue les types de données et les statistiques sommaires.

# Affichage des premières lignes de l'ensemble de données journalier
head(day_data)

## # A tibble: 6 × 16
##   instant dteday     season    yr  mnth holiday weekday workingday weathersit
##     <dbl> <date>      <dbl> <dbl> <dbl>   <dbl>   <dbl>      <dbl>      <dbl>
## 1       1 2011-01-01      1     0     1       0       6          0          2
## 2       2 2011-01-02      1     0     1       0       0          0          2
## 3       3 2011-01-03      1     0     1       0       1          1          1
## 4       4 2011-01-04      1     0     1       0       2          1          1
## 5       5 2011-01-05      1     0     1       0       3          1          1
## 6       6 2011-01-06      1     0     1       0       4          1          1
## # ℹ 7 more variables: temp <dbl>, atemp <dbl>, hum <dbl>, windspeed <dbl>,
## #   casual <dbl>, registered <dbl>, cnt <dbl>

# Résumé de l'ensemble de données journalier
summary(day_data)

##     instant          dteday               season            yr        
##  Min.   :  1.0   Min.   :2011-01-01   Min.   :1.000   Min.   :0.0000  
##  1st Qu.:183.5   1st Qu.:2011-07-02   1st Qu.:2.000   1st Qu.:0.0000  
##  Median :366.0   Median :2012-01-01   Median :3.000   Median :1.0000  
##  Mean   :366.0   Mean   :2012-01-01   Mean   :2.497   Mean   :0.5007  
##  3rd Qu.:548.5   3rd Qu.:2012-07-01   3rd Qu.:3.000   3rd Qu.:1.0000  
##  Max.   :731.0   Max.   :2012-12-31   Max.   :4.000   Max.   :1.0000  
##       mnth          holiday           weekday        workingday   
##  Min.   : 1.00   Min.   :0.00000   Min.   :0.000   Min.   :0.000  
##  1st Qu.: 4.00   1st Qu.:0.00000   1st Qu.:1.000   1st Qu.:0.000  
##  Median : 7.00   Median :0.00000   Median :3.000   Median :1.000  
##  Mean   : 6.52   Mean   :0.02873   Mean   :2.997   Mean   :0.684  
##  3rd Qu.:10.00   3rd Qu.:0.00000   3rd Qu.:5.000   3rd Qu.:1.000  
##  Max.   :12.00   Max.   :1.00000   Max.   :6.000   Max.   :1.000  
##    weathersit         temp             atemp              hum        
##  Min.   :1.000   Min.   :0.05913   Min.   :0.07907   Min.   :0.0000  
##  1st Qu.:1.000   1st Qu.:0.33708   1st Qu.:0.33784   1st Qu.:0.5200  
##  Median :1.000   Median :0.49833   Median :0.48673   Median :0.6267  
##  Mean   :1.395   Mean   :0.49538   Mean   :0.47435   Mean   :0.6279  
##  3rd Qu.:2.000   3rd Qu.:0.65542   3rd Qu.:0.60860   3rd Qu.:0.7302  
##  Max.   :3.000   Max.   :0.86167   Max.   :0.84090   Max.   :0.9725  
##    windspeed           casual         registered        cnt      
##  Min.   :0.02239   Min.   :   2.0   Min.   :  20   Min.   :  22  
##  1st Qu.:0.13495   1st Qu.: 315.5   1st Qu.:2497   1st Qu.:3152  
##  Median :0.18097   Median : 713.0   Median :3662   Median :4548  
##  Mean   :0.19049   Mean   : 848.2   Mean   :3656   Mean   :4504  
##  3rd Qu.:0.23321   3rd Qu.:1096.0   3rd Qu.:4776   3rd Qu.:5956  
##  Max.   :0.50746   Max.   :3410.0   Max.   :6946   Max.   :8714

# Structure de l'ensemble de données journalier
str(day_data)

## spc_tbl_ [731 × 16] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
##  $ instant   : num [1:731] 1 2 3 4 5 6 7 8 9 10 ...
##  $ dteday    : Date[1:731], format: "2011-01-01" "2011-01-02" ...
##  $ season    : num [1:731] 1 1 1 1 1 1 1 1 1 1 ...
##  $ yr        : num [1:731] 0 0 0 0 0 0 0 0 0 0 ...
##  $ mnth      : num [1:731] 1 1 1 1 1 1 1 1 1 1 ...
##  $ holiday   : num [1:731] 0 0 0 0 0 0 0 0 0 0 ...
##  $ weekday   : num [1:731] 6 0 1 2 3 4 5 6 0 1 ...
##  $ workingday: num [1:731] 0 0 1 1 1 1 1 0 0 1 ...
##  $ weathersit: num [1:731] 2 2 1 1 1 1 2 2 1 1 ...
##  $ temp      : num [1:731] 0.344 0.363 0.196 0.2 0.227 ...
##  $ atemp     : num [1:731] 0.364 0.354 0.189 0.212 0.229 ...
##  $ hum       : num [1:731] 0.806 0.696 0.437 0.59 0.437 ...
##  $ windspeed : num [1:731] 0.16 0.249 0.248 0.16 0.187 ...
##  $ casual    : num [1:731] 331 131 120 108 82 88 148 68 54 41 ...
##  $ registered: num [1:731] 654 670 1229 1454 1518 ...
##  $ cnt       : num [1:731] 985 801 1349 1562 1600 ...
##  - attr(*, "spec")=
##   .. cols(
##   ..   instant = col_double(),
##   ..   dteday = col_date(format = ""),
##   ..   season = col_double(),
##   ..   yr = col_double(),
##   ..   mnth = col_double(),
##   ..   holiday = col_double(),
##   ..   weekday = col_double(),
##   ..   workingday = col_double(),
##   ..   weathersit = col_double(),
##   ..   temp = col_double(),
##   ..   atemp = col_double(),
##   ..   hum = col_double(),
##   ..   windspeed = col_double(),
##   ..   casual = col_double(),
##   ..   registered = col_double(),
##   ..   cnt = col_double()
##   .. )
##  - attr(*, "problems")=<externalptr>

## # A tibble: 6 × 17
##   instant dteday     season    yr  mnth    hr holiday weekday workingday
##     <dbl> <date>      <dbl> <dbl> <dbl> <dbl>   <dbl>   <dbl>      <dbl>
## 1       1 2011-01-01      1     0     1     0       0       6          0
## 2       2 2011-01-01      1     0     1     1       0       6          0
## 3       3 2011-01-01      1     0     1     2       0       6          0
## 4       4 2011-01-01      1     0     1     3       0       6          0
## 5       5 2011-01-01      1     0     1     4       0       6          0
## 6       6 2011-01-01      1     0     1     5       0       6          0
## # ℹ 8 more variables: weathersit <dbl>, temp <dbl>, atemp <dbl>, hum <dbl>,
## #   windspeed <dbl>, casual <dbl>, registered <dbl>, cnt <dbl>

##     instant          dteday               season            yr        
##  Min.   :    1   Min.   :2011-01-01   Min.   :1.000   Min.   :0.0000  
##  1st Qu.: 4346   1st Qu.:2011-07-04   1st Qu.:2.000   1st Qu.:0.0000  
##  Median : 8690   Median :2012-01-02   Median :3.000   Median :1.0000  
##  Mean   : 8690   Mean   :2012-01-02   Mean   :2.502   Mean   :0.5026  
##  3rd Qu.:13034   3rd Qu.:2012-07-02   3rd Qu.:3.000   3rd Qu.:1.0000  
##  Max.   :17379   Max.   :2012-12-31   Max.   :4.000   Max.   :1.0000  
##       mnth              hr           holiday           weekday     
##  Min.   : 1.000   Min.   : 0.00   Min.   :0.00000   Min.   :0.000  
##  1st Qu.: 4.000   1st Qu.: 6.00   1st Qu.:0.00000   1st Qu.:1.000  
##  Median : 7.000   Median :12.00   Median :0.00000   Median :3.000  
##  Mean   : 6.538   Mean   :11.55   Mean   :0.02877   Mean   :3.004  
##  3rd Qu.:10.000   3rd Qu.:18.00   3rd Qu.:0.00000   3rd Qu.:5.000  
##  Max.   :12.000   Max.   :23.00   Max.   :1.00000   Max.   :6.000  
##    workingday       weathersit         temp           atemp       
##  Min.   :0.0000   Min.   :1.000   Min.   :0.020   Min.   :0.0000  
##  1st Qu.:0.0000   1st Qu.:1.000   1st Qu.:0.340   1st Qu.:0.3333  
##  Median :1.0000   Median :1.000   Median :0.500   Median :0.4848  
##  Mean   :0.6827   Mean   :1.425   Mean   :0.497   Mean   :0.4758  
##  3rd Qu.:1.0000   3rd Qu.:2.000   3rd Qu.:0.660   3rd Qu.:0.6212  
##  Max.   :1.0000   Max.   :4.000   Max.   :1.000   Max.   :1.0000  
##       hum           windspeed          casual         registered   
##  Min.   :0.0000   Min.   :0.0000   Min.   :  0.00   Min.   :  0.0  
##  1st Qu.:0.4800   1st Qu.:0.1045   1st Qu.:  4.00   1st Qu.: 34.0  
##  Median :0.6300   Median :0.1940   Median : 17.00   Median :115.0  
##  Mean   :0.6272   Mean   :0.1901   Mean   : 35.68   Mean   :153.8  
##  3rd Qu.:0.7800   3rd Qu.:0.2537   3rd Qu.: 48.00   3rd Qu.:220.0  
##  Max.   :1.0000   Max.   :0.8507   Max.   :367.00   Max.   :886.0  
##       cnt       
##  Min.   :  1.0  
##  1st Qu.: 40.0  
##  Median :142.0  
##  Mean   :189.5  
##  3rd Qu.:281.0  
##  Max.   :977.0

## spc_tbl_ [17,379 × 17] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
##  $ instant   : num [1:17379] 1 2 3 4 5 6 7 8 9 10 ...
##  $ dteday    : Date[1:17379], format: "2011-01-01" "2011-01-01" ...
##  $ season    : num [1:17379] 1 1 1 1 1 1 1 1 1 1 ...
##  $ yr        : num [1:17379] 0 0 0 0 0 0 0 0 0 0 ...
##  $ mnth      : num [1:17379] 1 1 1 1 1 1 1 1 1 1 ...
##  $ hr        : num [1:17379] 0 1 2 3 4 5 6 7 8 9 ...
##  $ holiday   : num [1:17379] 0 0 0 0 0 0 0 0 0 0 ...
##  $ weekday   : num [1:17379] 6 6 6 6 6 6 6 6 6 6 ...
##  $ workingday: num [1:17379] 0 0 0 0 0 0 0 0 0 0 ...
##  $ weathersit: num [1:17379] 1 1 1 1 1 2 1 1 1 1 ...
##  $ temp      : num [1:17379] 0.24 0.22 0.22 0.24 0.24 0.24 0.22 0.2 0.24 0.32 ...
##  $ atemp     : num [1:17379] 0.288 0.273 0.273 0.288 0.288 ...
##  $ hum       : num [1:17379] 0.81 0.8 0.8 0.75 0.75 0.75 0.8 0.86 0.75 0.76 ...
##  $ windspeed : num [1:17379] 0 0 0 0 0 0.0896 0 0 0 0 ...
##  $ casual    : num [1:17379] 3 8 5 3 0 0 2 1 1 8 ...
##  $ registered: num [1:17379] 13 32 27 10 1 1 0 2 7 6 ...
##  $ cnt       : num [1:17379] 16 40 32 13 1 1 2 3 8 14 ...
##  - attr(*, "spec")=
##   .. cols(
##   ..   instant = col_double(),
##   ..   dteday = col_date(format = ""),
##   ..   season = col_double(),
##   ..   yr = col_double(),
##   ..   mnth = col_double(),
##   ..   hr = col_double(),
##   ..   holiday = col_double(),
##   ..   weekday = col_double(),
##   ..   workingday = col_double(),
##   ..   weathersit = col_double(),
##   ..   temp = col_double(),
##   ..   atemp = col_double(),
##   ..   hum = col_double(),
##   ..   windspeed = col_double(),
##   ..   casual = col_double(),
##   ..   registered = col_double(),
##   ..   cnt = col_double()
##   .. )
##  - attr(*, "problems")=<externalptr>

sum(is.na(hour_data))

## [1] 0

sum(is.na(day_data))

## [1] 0

sum(duplicated(hour_data))

## [1] 0

sum(duplicated(day_data))

## [1] 0

Nous allons maintenant examiner nos données pour répondre aux questions spécifiques posées.

##Changements de Température Selon les Saisons #Comment les températures changent-elles selon les saisons ? Quelles sont les températures moyennes et médianes ?

Pour répondre à cette question, nous analyserons la température (temp) par saison et calculerons les moyennes et médianes.

library(dplyr)

# Calcul de la température moyenne et médiane par saison
temperature_stats <- day_data %>%
  group_by(season) %>%
  summarise(mean_temp = mean(temp, na.rm = TRUE),
            median_temp = median(temp, na.rm = TRUE))

# Affichage des résultats
print(temperature_stats)

## # A tibble: 4 × 3
##   season mean_temp median_temp
##    <dbl>     <dbl>       <dbl>
## 1      1     0.298       0.286
## 2      2     0.544       0.562
## 3      3     0.706       0.715
## 4      4     0.423       0.409

Cette partie du code permet d’obtenir une vue d’ensemble des températures par saison, ce qui est crucial pour comprendre comment les conditions météorologiques influencent l’utilisation des vélos.

##Y a-t-il une corrélation entre la température (temp/atemp) et le nombre total de locations de vélos ?

Nous examinerons la corrélation entre la température, la température ressentie et le nombre total de locations.

# Charger la bibliothèque corrplot si ce n'est pas déjà fait
# Si vous ne l'avez pas déjà installée, exécutez install.packages("corrplot")
library(corrplot)

## corrplot 0.92 loaded

# Calculer la température moyenne entre temp et atemp
day_data$mean_temp_atemp <- rowMeans(day_data[, c("temp", "atemp")])

# Calculer la matrice de corrélation
correlation_matrix <- cor(day_data[, c("temp", "atemp", "mean_temp_atemp", "cnt")], use = "complete.obs")

# Visualiser la matrice de corrélation
corrplot(correlation_matrix, method = "circle")

#Quelles sont les températures moyennes, l’humidité, la vitesse du vent et les locations totales par mois ?

Nous calculerons ces moyennes pour chaque mois pour identifier des tendances saisonnières.

Les valeurs de température, d’humidité et de vitesse du vent sont normalisées. La température moyenne est la plus élevée en juillet (0.755) et la plus basse en décembre (0.324). L’humidité est relativement stable tout au long de l’année avec une légère augmentation pendant les mois d’été. La vitesse du vent moyenne est généralement la plus faible en juillet (0.166) et la plus élevée en avril (0.234). Le nombre de locations de vélos est le plus élevé en juin (5772.37) et le plus bas en janvier (2176.34).

# Calcul des moyennes mensuelles
monthly_averages <- day_data %>%
  group_by(mnth) %>%
  summarise(
    mean_temp = mean(temp, na.rm = TRUE),
    mean_humidity = mean(hum, na.rm = TRUE),
    mean_windspeed = mean(windspeed, na.rm = TRUE),
    mean_rentals = mean(cnt, na.rm = TRUE)
  )

# Affichage des moyennes mensuelles
print(monthly_averages)

## # A tibble: 12 × 5
##     mnth mean_temp mean_humidity mean_windspeed mean_rentals
##    <dbl>     <dbl>         <dbl>          <dbl>        <dbl>
##  1     1     0.236         0.586          0.206        2176.
##  2     2     0.299         0.567          0.216        2655.
##  3     3     0.391         0.588          0.223        3692.
##  4     4     0.470         0.588          0.234        4485.
##  5     5     0.595         0.689          0.183        5350.
##  6     6     0.684         0.576          0.185        5772.
##  7     7     0.755         0.598          0.166        5564.
##  8     8     0.709         0.638          0.173        5664.
##  9     9     0.616         0.715          0.166        5767.
## 10    10     0.485         0.694          0.175        5199.
## 11    11     0.369         0.625          0.184        4247.
## 12    12     0.324         0.666          0.177        3404.

#La température est-elle associée aux locations de vélos (enregistrés vs occasionnels) ?

Nous étudierons comment la température affecte différemment les utilisateurs enregistrés et occasionnels.

# Association de la température avec les locations de vélos
rentals_by_temp <- day_data %>%
  group_by(temp) %>%
  summarise(mean_casual = mean(casual),
            mean_registered = mean(registered))

# Visualisation
ggplot(rentals_by_temp, aes(x = temp)) +
  geom_line(aes(y = mean_casual, color = "Occasionnel")) +
  geom_line(aes(y = mean_registered, color = "Enregistré")) +
  labs(title = "Locations de Vélos par Température",
       x = "Température Normalisée",
       y = "Locations Moyennes")

# Calcul de la corrélation entre température et locations pour les utilisateurs enregistrés et occasionnels
correlation_temp_users <- day_data %>%
  select(temp, casual, registered) %>%
  cor(use = "complete.obs")

# Affichage de la matrice de corrélation
print(correlation_temp_users)

##                 temp    casual registered
## temp       1.0000000 0.5432847  0.5400120
## casual     0.5432847 1.0000000  0.3952825
## registered 0.5400120 0.3952825  1.0000000

Calculer la corrélation entre la température et les locations pour les utilisateurs enregistrés et occasionnels.

Résumer comment les locations moyennes varient avec la température en regroupant les températures par intervalles puis en calculant le nombre moyen de locations pour chaque groupe. ????????????????????je ne sais pas

# Conversion de 'dteday' en format de date
day_data$dteday <- as.Date(day_data$dteday, format="%Y-%m-%d")

# Tracé de 'cnt' vs 'dteday'
ggplot(day_data, aes(x=dteday, y=cnt)) + 
  geom_line() + 
  labs(title="Nombre de locations de vélos par jour", 
       x="Date", 
       y="Nombre total de locations de vélos") +
  theme(axis.text.x = element_text(angle=90, hjust=1))

Il y a une saison annuelle, on voit une tendance ( mois plus chauds implique hausse des locations et une baisse pendant les mois plus froids). La variance est plus au moins constante.

# Conversion de 'dteday' en format de date
day_data$dteday <- as.Date(day_data$dteday, format="%Y-%m-%d")

# Identifier les valeurs aberrantes, par exemple en utilisant la méthode des écarts interquartiles
IQR_values <- IQR(day_data$cnt)
quantiles <- quantile(day_data$cnt, probs=c(.25, .75))
cap <- quantiles[2] + 1.5 * IQR_values
floor <- quantiles[1] - 1.5 * IQR_values

# Filtrer les valeurs aberrantes
day_data <- day_data %>%
  filter(cnt >= floor & cnt <= cap)

# Lisser la série temporelle en utilisant une moyenne mobile avec la fonction 'filter' du package 'stats'
day_data$cnttimeseries <- ts(day_data$cnt, frequency = 12)

hw_model <- HoltWinters(day_data$cnttimeseries)

plot(hw_model)

# Extraire les valeurs prévues à partir du modèle HoltWinters
fitted_values <- hw_model$fitted[, "xhat"]

# Ajouter la fréquence appropriée à la série lissée (remplacez 12 par la fréquence appropriée)
smoothed_ts <- ts(fitted_values, frequency = 12)

plot(smoothed_ts)

result <- adf.test(smoothed_ts, alternative = "stationary")

# Afficher les résultats du test
print(result)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  smoothed_ts
## Dickey-Fuller = -0.7527, Lag order = 8, p-value = 0.9658
## alternative hypothesis: stationary

La série n’est pas stationnaire (variance constante mais moyenne variable) le test de Dicker Fuller confirme cela avec une p_value = 0.96 (hyp nulle : non stationnaire), on distingue une saison (qui se répète deux fois).

# Appliquer la différenciation saisonnière
diff_series <- diff(smoothed_ts, lag=12)

# Tracer la série différenciée
plot(diff_series, main="Série Temporelle Différenciée", xlab="Temps", ylab="Différences")

pacf(diff_series,lag.max=100, main="PACF")

acf(diff_series,lag.max= 100, main="ACF")

On remarque le PACF décroit exponentiellement, et sur le ACF on lit : MA(9) On remarque le ACF décroit exponentiellement, et sur le PACF on lit : AR(2)

Concernant les saisons, sur le ACF ça décroit et sur le PACF on lit P = 4 sur le PACF ça décroit et sur le ACF on lit Q=4

#MA(9) ou AR(2)
#ARIMA(2,0,0)(3,1,0)h=12 ***
#ARIMA(0,0,9)(4,1,0)h=12 ****
#ARIMA(2,0,0)(0,1,2)h=12 *
#ARIMA(0,0,9)(0,1,2)h=12 **

model_m1 <- arima(smoothed_ts, order = c(2, 0, 0), seasonal = list(order = c(3, 1, 0), period = 12))

residuals_m1 <- residuals(model_m1)
plot(residuals_m1)

#mean=0 var=cte 
par(mfrow = c(2, 1))
acf(residuals_m1, main = "ACF of Residuals")
pacf(residuals_m1, main = "PACF of Residuals")

#no significant peak on the ACF and PACF

Box.test(residuals_m1, type = "Ljung-Box")

## 
##  Box-Ljung test
## 
## data:  residuals_m1
## X-squared = 0.057693, df = 1, p-value = 0.8102

#p-value= 0.8102
aic_value_m1 <- AIC(model_m1)
cat("AIC du modèle m1:", aic_value_m1, "\n")

## AIC du modèle m1: 10179.9

#AIC=10179.9 

model_m2 <- arima(smoothed_ts, order = c(0, 0, 9), seasonal = list(order = c(4, 1, 0), period = 12))

residuals_m2 <- residuals(model_m2)
plot(residuals_m2)

#mean=0 var=cte 
par(mfrow = c(2, 1))
acf(residuals_m2, main = "ACF of Residuals")
pacf(residuals_m2, main = "PACF of Residuals")

#no significant peak on the ACF and PACF

Box.test(residuals_m2, type = "Ljung-Box")

## 
##  Box-Ljung test
## 
## data:  residuals_m2
## X-squared = 0.064682, df = 1, p-value = 0.7992

#p-value= 0.7992
aic_value_m2 <- AIC(model_m2)
cat("AIC du modèle m2:", aic_value_m2, "\n")

## AIC du modèle m2: 10213.16

#AIC=10213.16 

model_m3 <- arima(smoothed_ts, order = c(2, 0, 0), seasonal = list(order = c(0, 1, 2), period = 12))

residuals_m3 <- residuals(model_m3)
plot(residuals_m3)

#mean=0 var=cte 
par(mfrow = c(2, 1))
acf(residuals_m3, main = "ACF of Residuals")
pacf(residuals_m3, main = "PACF of Residuals")

#no significant peak on the ACF and PACF

Box.test(residuals_m3, type = "Ljung-Box")

## 
##  Box-Ljung test
## 
## data:  residuals_m3
## X-squared = 0.126, df = 1, p-value = 0.7226

#p-value= 0.7226
aic_value_m3 <- AIC(model_m3)
cat("AIC du modèle m3:", aic_value_m3, "\n")

## AIC du modèle m3: 10131.05

#AIC=10131.05  

model_m4 <- arima(smoothed_ts, order = c(0, 0, 9), seasonal = list(order = c(0, 1, 2), period = 12))

residuals_m4 <- residuals(model_m4)
plot(residuals_m4)

#mean=0 var=cte 
par(mfrow = c(2, 1))
acf(residuals_m4, main = "ACF of Residuals")
pacf(residuals_m4, main = "PACF of Residuals")

#no significant peak on the ACF and PACF

Box.test(residuals_m4, type = "Ljung-Box")

## 
##  Box-Ljung test
## 
## data:  residuals_m4
## X-squared = 0.051688, df = 1, p-value = 0.8202

#p-value= 0.8202
aic_value_m4 <- AIC(model_m4)
cat("AIC du modèle m4:", aic_value_m4, "\n")

## AIC du modèle m4: 10206.35

#AIC=10206.35 

Le meilleur modèle a l’air d’être le modèle m2 : ARIMA(0,0,9)(4,1,0)h=12

library(forecast)

# Ajuster un modèle ARIMA automatiquement sur les données stationnaires
# Supposons que 'diff_log_cnt_clean' est votre série temporelle prétraitée
model <- auto.arima(smoothed_ts)


# Afficher le résumé du modèle ajusté
summary(model)

## Series: smoothed_ts 
## ARIMA(5,1,1)(2,0,0)[12] 
## 
## Coefficients:
##           ar1      ar2      ar3      ar4      ar5     ma1    sar1    sar2
##       -0.5396  -0.0782  -0.1907  -0.2164  -0.1897  0.6284  0.2808  0.3158
## s.e.   0.2044   0.0463   0.0472   0.0470   0.0402  0.2054  0.0372  0.0370
## 
## sigma^2 = 104406:  log likelihood = -5165.91
## AIC=10349.81   AICc=10350.07   BIC=10391
## 
## Training set error measures:
##                     ME     RMSE      MAE        MPE     MAPE      MASE
## Training set -3.822976 321.0912 235.8936 -0.6894523 6.372798 0.4056763
##                     ACF1
## Training set 0.002163913

##Désaisonnaliser les données :

#Désaisonnaliser les données 
diff_cnt <- diff(day_data$cnt, lag=12)
plot(diff_cnt)

pacf(diff_cnt, main="PACF")

L’ACF présente une décroissance rapide après le premier décalage, ce qui suggère un terme MA. Si ce décalage est le seul significatif, cela pourrait indiquer un modèle MA(1).

Le PACF montre un pic significatif au premier décalage et ensuite il se stabilise, ce qui est typique d’un modèle AR(1).

acf(diff_cnt, main="ACF")

#AR(1)
#MA(3)
#ARIMA(1,0,0)
#ARIMA(0,0,3)

model_m1_deseasonal <- arima(diff_cnt, order = c(1, 0, 0))

residuals_m1_deseasonal <- residuals(model_m1_deseasonal)
plot(residuals_m1_deseasonal)

#mean=0 var=cte 
par(mfrow = c(2, 1))
acf(residuals_m1_deseasonal, main = "ACF of Residuals")
pacf(residuals_m1_deseasonal, main = "PACF of Residuals")

#no significant peak on the ACF and PACF

Box.test(residuals_m1_deseasonal, type = "Ljung-Box")

## 
##  Box-Ljung test
## 
## data:  residuals_m1_deseasonal
## X-squared = 0.28091, df = 1, p-value = 0.5961

#p-value= 0.5961
aic_value_m1_deseasonal<- AIC(model_m1_deseasonal)
cat("AIC du modèle m4:", aic_value_m1_deseasonal, "\n")

## AIC du modèle m4: 12343.73

#AIC=12343.73  

model_m2_deseasonal <- arima(diff_cnt, order = c(0, 0, 3))

residuals_m2_deseasonal <- residuals(model_m2_deseasonal)
plot(residuals_m2_deseasonal)

#mean=0 var=cte 
par(mfrow = c(2, 1))
acf(residuals_m2_deseasonal, main = "ACF of Residuals")
pacf(residuals_m2_deseasonal, main = "PACF of Residuals")

#no significant peak on the ACF and PACF

Box.test(residuals_m2_deseasonal, type = "Ljung-Box")

## 
##  Box-Ljung test
## 
## data:  residuals_m2_deseasonal
## X-squared = 0.0016148, df = 1, p-value = 0.9679

#p-value= 0.9679
aic_value_m2_deseasonal<- AIC(model_m2_deseasonal)
cat("AIC du modèle m2:", aic_value_m2_deseasonal, "\n")

## AIC du modèle m2: 12346.61

#AIC=12346.61    

Le meilleur modèle en prenant en compte the Parsimony principle est le modèle m1 : AR(1)

# Installer le package forecast si nécessaire

# Charger le package forecast
library(forecast)

# Ajuster un modèle ARIMA automatiquement sur les données stationnaires
fit_arima <- auto.arima(diff_cnt)

# Afficher le résumé du modèle ajusté
summary(fit_arima)

## Series: diff_cnt 
## ARIMA(0,1,3) 
## 
## Coefficients:
##           ma1      ma2      ma3
##       -0.5378  -0.3000  -0.1367
## s.e.   0.0367   0.0414   0.0361
## 
## sigma^2 = 1664577:  log likelihood = -6161.37
## AIC=12330.74   AICc=12330.8   BIC=12349.05
## 
## Training set error measures:
##                     ME     RMSE      MAE      MPE     MAPE      MASE
## Training set -44.02235 1286.591 924.9182 97.16658 163.5412 0.8277705
##                     ACF1
## Training set 0.001903797

residuals_adjusted_deseasonal <- residuals(fit_arima)
plot(residuals_adjusted_deseasonal)

#mean=0 var=cte 
par(mfrow = c(2, 1))
acf(residuals_adjusted_deseasonal, main = "ACF of Residuals")
pacf(residuals_adjusted_deseasonal, main = "PACF of Residuals")

#no significant peak on the ACF and PACF

Box.test(residuals_adjusted_deseasonal, type = "Ljung-Box")

## 
##  Box-Ljung test
## 
## data:  residuals_adjusted_deseasonal
## X-squared = 0.0026169, df = 1, p-value = 0.9592

aic_value_adjusted<- AIC(fit_arima)
cat("AIC du modèle auto arima:", aic_value_adjusted, "\n")

## AIC du modèle auto arima: 12330.74

Le modèle choisi plus haut, AR(1) donne un AIC légèrement plus élevé, de plus, il est plus simple. et d’après le principe de parsimony, il est préférable de choisir le plus simple.

library(forecast)
forecasts <- forecast(model_m1_deseasonal, h=25)
plot(forecasts)

cat(length(day_data$cnt))

## 731

# Diviser les données en ensembles d'entraînement et de test
train_data <- window(day_data$cnt, start = 1, end = 700)
test_data <- window(day_data$cnt, start = 701)

train_data_ts <- ts(train_data)
test_data_ts <- ts(test_data)

ts.plot(train_data_ts)

diff_train <- diff(train_data_ts, frequency=12)
plot(diff_train)

pacf(diff_train, main="PACF")

acf(diff_train, main="ACF")

On remarque ACF qui décroit donc AR(5) PACF décroît donc MA(3)

auto_model <- auto.arima(train_data)
summary(auto_model)

## Series: train_data 
## ARIMA(3,1,1) 
## 
## Coefficients:
##          ar1      ar2      ar3      ma1
##       0.3160  -0.0401  -0.0916  -0.8718
## s.e.  0.0437   0.0410   0.0409   0.0237
## 
## sigma^2 = 831740:  log likelihood = -5754.54
## AIC=11519.08   AICc=11519.16   BIC=11541.82
## 
## Training set error measures:
##                    ME     RMSE      MAE       MPE     MAPE     MASE
## Training set 32.40763 908.7347 637.0445 -44.80348 59.12826 0.876924
##                      ACF1
## Training set -0.006730811

# Ajuster un modèle ARIMA manuellement sur l'ensemble d'entraînement
manual_model <- arima(train_data, order = c(0, 0, 3), seasonal = list(order = c(4, 1, 0), period = 12))

# Ajuster un modèle ARIMA automatiquement sur l'ensemble d'entraînement
auto_model <- auto.arima(train_data)

# Prévoir les prochaines 25 observations avec les deux modèles
manual_forecast <- forecast(manual_model, h = 25)
auto_forecast <- forecast(auto_model, h = 25)

# Plot des prévisions et des données réelles
plot(test_data, col = "blue", ylim = range(test_data, manual_forecast$mean, auto_forecast$mean),
     main = "ARIMA Forecast Comparison", xlab = "Time", ylab = "Value")
lines(manual_forecast$mean, col = "red", lty = 2, lwd = 2)
lines(auto_forecast$mean, col = "green", lty = 2, lwd = 2)

# Ajouter une légende
legend("topright", legend = c("Observed", "Manual Forecast", "Auto Forecast"),
       col = c("blue", "red", "green"), lty = c(1, 2, 2), lwd = 2)